/*
    Copyright 2008-2011 Patrik Jonsson, patrik@familjenjonsson.org,
    Chris Hayward

    This file is part of Sunrise.

    Sunrise is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 3 of the License, or
    (at your option) any later version.

    Sunrise is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with Sunrise.  If not, see <http://www.gnu.org/licenses/>.

*/

/// \file
/// Contains stellarmodel, a class for reading a stellar model FITS file as
/// generated by convertstellarmodel. \ingroup sfrhist

// $Id$

#include "bhmodel.h"
#include "preferences.h"
#include <memory> 
#include <fstream> 
#include "CCfits/CCfits"
#include "blitz-fits.h"
#include "fits-utilities.h"
#include "misc.h"
#include "interpolatort.h"
#include "constants.h"
#include "mcrx-debug.h"
#include "vecops.h"

using namespace std;
using namespace CCfits;
using namespace blitz;

mcrx::bhmodel::bhmodel(Preferences& p) : prefs(p)

{
  // The first thing we should do is to open the FITS file containing
  // the black hole model specified in keyword "bhmodelfile"
  auto_ptr<FITS> input_file;

  if (prefs.defined("bhmodelfile")) {
    file_name = 
      word_expand(prefs.getValue("bhmodelfile", string ()))[0];
    try {
      input_file.reset(new FITS (file_name, Read));//, "SED", true));
    }
    catch (...) {
      cerr << "Error: Can't open black hole model file: "
           << file_name << endl;
      exit (1); 
    }
  }
  else {
    cerr << "Error: No keyword \"bhmodelfile\" defined" <<endl;
    exit (1);
  }
  cout << "Reading black hole model file: " << file_name <<endl;

  // Check that the file is actually a black hole model file
  string file_type;
  input_file->pHDU().readKey ("FILETYPE", file_type);
  if (file_type != "BHMODEL") {
    cerr << "Error: File " << file_type << " is not a black hole model file"
         <<endl;
    exit (1);
  } 

  ExtHDU & parameters = open_HDU(*input_file,"BHMODEL");
  parameters.readKey("logflux",  log_flux);

  // Read the SED
  ExtHDU & spectrum = open_HDU(*input_file, "SED");
  read (spectrum, sed);
  

  // Read luminosities.
  ExtHDU & l = open_HDU (*input_file, "AXES");
  l.column("log_L_bol").read(log_L_bol, 1, n_lum () );
  l.column("lambda").read(lambdas, 1, n_lambda () );

  units_["luminosity"] = l.column("log_L_bol").unit();
  units_["wavelength"] = l.column("lambda").unit();

  // Check that the units are what we expect
  if ((units().get("luminosity") != "W") || 
      (units().get("wavelength") != "m")) {
    cerr << "WARNING: Nonstandard units in black hole model file\n\t" 
	 << units().get("luminosity") << '\t' << units().get("wavelengths")
	 << endl;
  }
  // Check if FITS file contains log flux; if so, convert to flux.
  if (log_flux)
     sed = pow (10., sed);

  // Update to reflect correction.
  log_flux = false;

  // check energy of BH seds
  alambdas.reference(reference_vector(lambdas));

  // make sure the SED extends over the requested range, if not pad
  // with zero on the short end and a 1/lambda^2 decline on the long
  // end.  All these indices are getting very messy and difficult to
  // keep track of, so we should clean this up.

  // Note that we should NEVER truncate, because it will be resampled
  // anyway. We just want to make sure it doesn't fail if the
  // wavelength range asked for is larger.

  // We choose the points extend ACROSS the extended range in order to
  // avoid problems with resampling.
  bool extend_front=false;
  bool extend_back=false;
  if(p.getValue("min_wavelength", T_float()) < lambdas.front()) {
    lambdas.insert(lambdas.begin(), 
		   p.getValue("min_wavelength", T_float())*0.999);
    extend_front=true;
  }

  if(p.getValue("max_wavelength", T_float()) > lambdas.back()) {
    lambdas.push_back(p.getValue("max_wavelength", T_float())*1.001);
    extend_back=true;
  }

  const int new_lambda_len = lambdas.size() + (extend_front?1:0) + (extend_back?1:0);

  array_2 new_sed(n_lum(), new_lambda_len);
  new_sed(Range::all(), Range(extend_front?1:0, sed.extent(secondDim) + (extend_front?1:0) - 1 ) ) = 
    sed(Range::all(), Range::all());

  if(extend_front)
    new_sed(Range::all(), 0) = blitz::tiny(T_float());
  if(extend_back)
    // calculate end point as L_lambda ~ 1/lambda^4, assuming we are on
    // the Rayleigh-Jeans tail of the SED
    new_sed(Range::all(), new_lambda_len-1) =
      new_sed(Range::all(), new_lambda_len-2)*
      pow(lambdas[lambdas.size()-2]/lambdas.back(), 4);

  sed.reference(new_sed);

  alambdas.reference(reference_vector(lambdas));

  middle_log_L_bol = calculate_middle(log_L_bol);

  //convert to bolometric luminosity in W; must convert from (M_solar/yr)
  //to (kg/s) then use L = eff* mdotbh c^2
  /// \todo NEED TO TAKE UNITS INTO ACCOUNT HERE
  mdot2lum_ = prefs.getValue("bh_radiative_efficiency", T_float())*
    constants::c0*constants::c0*constants::Msun/constants::yr;
}



void
mcrx::bhmodel::write_fits_parameters (const std::string& output_file_name) const
{
  FITS output_file(output_file_name, Write);
  auto_ptr<FITS> input_file (new FITS (file_name, Read));

  ExtHDU& st=open_HDU(*input_file,"BHMODEL");
  output_file.copy(st);
}

mcrx::array_1 mcrx::bhmodel::get_SED(bh_particle_set& s, int index) const
{
  //get BH mass accretion rate
  const T_float mdotbh = s.mdotbh(index);
  const T_float L_bol = mdot2lum_*mdotbh;
  cout << "Black hole mdot = " << mdotbh << " L_bol = "<< L_bol << endl;

  //find L_bol index
  const int l = (L_bol == 0) ? 0 :
		find_if (middle_log_L_bol.begin(), middle_log_L_bol.end(),
			 bind2nd (greater<double> (), log10(L_bol))) -
		middle_log_L_bol.begin() - 1;
  assert (l >= 0);
  assert (l < n_lum());

  //return SED array for L_bol index l, renormalized to L_bol of BH
  DEBUG(1,cout << "BH SED " << l << " has bolometric luminosity " << pow(10.,log_L_bol.at (l)) << endl;);

  /// \todo THIS IS INCORRECT IF LOGFLUX==TRUE. GET RID OF LOGFLUX.
  assert(!log_flux);
  const T_float renorm = L_bol / pow(10.,log_L_bol[l]);
  return array_1(sed (l, Range::all () ) * renorm);
}


void mcrx::bhmodel::resample(const array_1& new_lambda)
{
  // check if we are upsampling or downsampling. If upsampling, we use
  // mcrx::resample with a logarithmic resampling, because that works
  // better. For downsampling, we use mcrx::subsample, which preserves
  // the integral.
  const bool upsample = 
    (log10(lambda()(lambda().size()-1))-log10(lambda()(0)))/(lambda().size()-1) >
    ((log10(new_lambda(new_lambda.size()-1))-log10(new_lambda(0)))/
     (new_lambda.size()-1));
  if (upsample) 
    cout << "Upsampling SEDs onto wavelength grid\n";
  else
    cout << "Downsampling SEDs onto wavelength grid\n";

  array_2 new_sed(sed.extent(firstDim),new_lambda.size());

  if (upsample)
    for(int i=0; i <sed.extent(firstDim); ++i) {
      DEBUG(1,cout << "bh sed " << i << '\t'	     \
	    << integrate_quantity(sed(i,Range::all()),	\
				  alambdas)		\
	    << '\t' << pow(10.,log_L_bol[i]) << endl;);
      new_sed(i,Range::all()) = mcrx::resample(sed(i,Range::all()),
					       alambdas, new_lambda, 
					       false, true);
    }
  else {
    for(int i=0; i <sed.extent(firstDim); ++i) {
      DEBUG(1,cout << "bh sed " << i << '\t'	     \
	    << integrate_quantity(sed(i,Range::all()),	\
				  alambdas)		\
	    << '\t' << pow(10.,log_L_bol[i]) << endl;);
      // use subsample to ensure we preserve luminosity
      new_sed(i,Range::all()) = mcrx::subsample(sed(i,Range::all()),
						alambdas, new_lambda);
    }

  }
  sed.reference(new_sed);
  lambdas.assign(new_lambda.begin(), new_lambda.end());
  alambdas.reference(reference_vector(lambdas));
}
